using System;

namespace UraniaLib
{
	/// <summary>
	/// Summary description for UraniaSunriseSunset.
	/// </summary>
	public class UraniaSunriseSunset
	{
		public UraniaSunriseSunset()
		{
			//
			// TODO: Add constructor logic here
			//
		}


		//
		//   This module contains a user defined spreadsheet function called
		//   sunevent(year, month, day, glong, glat, tz,event) that returns the
		//   time of day of sunrise, sunset, or the beginning and end one of
		//   three kinds of twilight
		//
		//   The method used here is adapted from Montenbruck and Pfleger's
		//   Astronomy on the Personal Computer, 3rd Ed, section 3.8
		//
		//   the arguments for the function are as follows...
		//   year, month, day - your date in zone time
		//   glong - your longitude in degrees, west negative
		//   glat - your latitude in degrees, north positive
		//   tz - your time zone in decimal hours, west or 'behind' Greenwich negative
		//   event - a code integer representing the event you want as follows
		//       positive integers mean rising events
		//       negative integers mean setting events
		//       1 = sunrise                 -1  = sunset
		//       2 = begin civil twilight    -2  = end civil twilight
		//       3 = begin nautical twi      -3  = end nautical twi
		//       4 = begin astro twi         -4  = end astro twi
		//
		//   the results are returned as a variant with either a time of day
		//   in zone time or a string reporting an 'event not occuring' condition
		//   event not occuring can be one of the following
		//       .....    always below horizon, so no rise or set
		//       *****    always above horizon, so no rise or set
		//       -----    the particular rise or set event does not occur on that day
		//
		//   The function will produce meaningful results at all latitudes
		//   but there will be a small range of latitudes around 67.43 degrees North or South
		//   when the function might indicate a sunrise very close to noon (or a sunset
		//   very soon after noon) where in fact the Sun is below the horizon all day
		//   this behaviour relates to the approximate Sun position formulas in use
		//
		//   As always, the sunrise / set times relate to an earth which is smooth
		//   and has no large obstructions on the horizon - you might get a close
		//   approximation to this at sea but rarely on land. Accuracy more than 1 min
		//   of time is not worth striving for - atmospheric refraction alone can
		//   alter observed rise times by minutes
		//
		//   The module also defines a series of 'named funtions' based on sunevent()
		//   as follows
		//   astrotwilightstarts(date, tz, glong, glat)
		//   nauticaltwilightstarts(date, tz, glong, glat)
		//   civiltwilightstarts(date, tz, glong, glat)
		//   sunrise(date, tz, glong, glat)
		//   sunset(date, tz, glong, glat)
		//   civiltwilightends(date, tz, glong, glat)
		//   nauticaltwilightends(date, tz, glong, glat)
		//   astrotwilightends(date, tz, glong, glat)
		//
		//   these functions  take a date in Excel date format and return times in
		//   Excel time format (ie a day fraction) if there is an event on the day
		//   If there isn't an event on the day, then you get #VALUE!
		//   These functions lend themselves to plotting sunrise and set times on
		//   charts - just multiply output in 'number' format by 24 to get the decimal
		//   hours for the event
		//
		
		public static int cSunrise = 1;
		public static int cBeginCivilTwilight = - 2;
		public static int cBeginNautTwilight = - 3;
		public static int cBeginAstroTwilight = - 4;
		public static int cSunset = - 1;
		public static int cEndCivilTwilight = 2;
		public static int cEndNautTwilight = 3;
		public static int cEndAstroTwilight = 4;
		
		
		private static double mjd(int year, int month, int day)
		{
			//
			//   takes the year, month and day as a Gregorian calendar date
			//   and returns the modified julian day number
			//
			double a;
			double b;
			
			if (month <= 2)
			{
				month = month + 12;
				year--;
			}
			
			b = (int)Math.Floor(year / 400.0) - (int)Math.Floor(year / 100.0) + (int)Math.Floor(year / 4.0);
			a = 365.0 * year - 679004.0;
			return a + b + (int)Math.Floor(30.6001 *(month + 1)) + day;
		}
		
		
		private static double frac(double x)
		{
			//
			//  returns the fractional part of x as used in minimoon and minisun
			//
			double a;
			a = x - (int)Math.Floor(x);
			return a;
		}
		
		private static double range(double x)
		{
			//
			//   returns an angle in degrees in the range 0 to 360
			//   used to condition the arguments for the Sun's orbit
			//   in function minisun below
			//
			double a;
			double b;
			b = x / 360;
			a = 360 *(b - (int)Math.Floor(b));
			if (a < 0)
			{
				a = a + 360;
			}
			return a;
		}
		
		private static string hrsmin(double t)
		{
			//
			//   takes a time as a decimal number of hours between 0 and 23.9999...
			//   and returns a string with the time in hhmm format
			//
			double hour;
			double min;
			
			hour = (int)Math.Floor(t);
			min = (int)Math.Floor((t - hour) * 60 + 0.5);
			return (hour.ToString("00") + min.ToString("00"));
			//return "0000";
		}
		
		
		private static double lmst(double mjd, double glong)
		{
			//
			//  Takes the mjd and the longitude (west negative) and then returns
			//  the local sidereal time in hours. Im using Meeus formula 11.4
			//  instead of messing about with UTo and so on
			//
			double lst;
			double t;
			double d;
			d = mjd - 51544.5;
			t = d / 36525.0;
			lst = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - t * t * t / 38710000.0);
			return lst / 15.0 + glong / 15.0;
		}
		
		private static void minisun(double t, ref double ra, ref double dec)
		{
			//
			//   takes t (julian centuries since J2000.0) and empty variables ra and dec
			//   sets ra and dec to the value of the Sun coordinates at t
			//
			//   positions claimed to be within 1 arc min by Montenbruck and Pfleger
			//
			double p2;
			double coseps;
			double sineps;
			double l;
			double m;
			double DL;
			double SL;
			double x;
			double y;
			double z;
			double rho;
			p2 = 6.283185307;
			coseps = 0.91748;
			sineps = 0.39778;
			
			m = p2 * frac(0.993133 + 99.997361 * t);
			DL = 6893.0 * System.Math.Sin(m) + 72.0 * System.Math.Sin(2 * m);
			l = p2 * frac(0.7859453 + m / p2 +(6191.2 * t + DL) / 1296000.0);
			SL = System.Math.Sin(l);
			x = System.Math.Cos(l);
			y = coseps * SL;
			z = sineps * SL;
			rho = System.Math.Sqrt(1 - z * z);
			dec = (360.0 / p2) * System.Math.Atan(z / rho);
			ra = (48.0 / p2) * System.Math.Atan(y /(x + rho));
			if (ra < 0)
			{
				ra = ra + 24;
			}
		}
		
		
		private static void quad(double ym, double yz, double yp, ref int nz, ref double z1, ref double z2, ref double xe, ref double ye)
		{
			//
			//  finds the parabola throuh the three points (-1,ym), (0,yz), (1, yp)
			//  and sets the coordinates of the max/min (if any) xe, ye
			//  the values of x where the parabola crosses zero (z1, z2)
			//  and the nz number of roots (0, 1 or 2) within the interval [-1, 1]
			//
			double a;
			double b;
			double c;
			double dis;
			double dx;
			
			nz = 0;
			a = 0.5 *(ym + yp) - yz;
			b = 0.5 *(yp - ym);
			c = yz;
			xe = - b /(2 * a);
			ye = (a * xe + b) * xe + c;
			dis = b * b - 4.0 * a * c;
			if (dis > 0)
			{
				dx = 0.5 * System.Math.Sqrt(dis) / System.Math.Abs(a);
				z1 = xe - dx;
				z2 = xe + dx;
				if (System.Math.Abs(z1) <= 1.0)
				{
					nz++;
				}
				if (System.Math.Abs(z2) <= 1.0)
				{
					nz++;
				}
				if (z1 < - 1.0)
				{
					z1 = z2;
				}
			}
		}
		
		private static double SinAltSun(double mjd0, double hour, double glong, double cglat, double sglat)
		{
			//
			//  this rather mickey mouse function takes a lot of
			//  arguments and then returns the sine of the altitude of
			//  the object labelled by iobj. iobj = 1 is moon, iobj = 2 is sun
			//
			double mjd;
			double t;
			double ra = 0;
			double dec = 0;
			double tau;
			double salt;
			double rads;
			rads = 0.0174532925;
			mjd = mjd0 + hour / 24.0;
			t = (mjd - 51544.5) / 36525.0;
			minisun(t, ref ra, ref dec);
			// hour angle of object
			tau = 15.0 *(lmst(mjd, glong) - ra);
			// sin(alt) of object using the conversion formulas
			salt = sglat * System.Math.Sin(rads * dec) + cglat * System.Math.Cos(rads * dec) * System.Math.Cos(rads * tau);
			return salt;
		}
		
		//
		//   Worksheet functions below....
		//
		
		public static string sunevent(int year, int month, int day, double tz, double glong, double glat, int EventType)
		{
			//
			//   This is the function that does most of the work
			//
			//			double sglong;
			double sglat;
			double cglat;
			double ddate;
			double ym;
			double yz;
			int above;
			double utrise = 0;
			double utset = 0;
			// 			int above;
			utrise = 0;
			// 			double utset;
			// 			int above;
			// 			double utrise;
			// 			double utset;
			double yp;
			int nz;
			int rise;
			int sett;
			int j;
			double hour;
			double z1;
			double z2;
			double rads;
			double xe;
			double ye;
			string AlwaysUp;
			string AlwaysDown;
			string OutString = "";
			string NoEvent;
			// 			string AlwaysDown;
			// 			string OutString;
			// 			string NoEvent;
			double[] sinho = new double[6];
			rads = 0.0174532925;
			AlwaysUp = "****";
			AlwaysDown = "....";
			NoEvent = "----";
			
			//
			//   Set up the array with the 4 values of sinho needed for the 4
			//   kinds of sun event
			//
			sinho[1] = System.Math.Sin(rads * - 0.833); //sunset upper limb simple refraction
			sinho[2] = System.Math.Sin(rads * - 6.0); //civil twi
			sinho[3] = System.Math.Sin(rads * - 12.0); //nautical twi
			sinho[4] = System.Math.Sin(rads * - 18.0); //astro twi
			sglat = System.Math.Sin(rads * glat);
			cglat = System.Math.Cos(rads * glat);
			ddate = mjd(year, month, day) - tz / 24.0;
			//
			//   main loop takes each value of sinho in turn and finds the rise/set
			//   events associated with that altitude of the Sun
			//
			j = System.Math.Abs(EventType);
			nz = 0;
			z1 = 0;
			z2 = 0;
			xe = 0;
			ye = 0;
			rise = 0;
			sett = 0;
			above = 0;
			hour = 1.0;
			ym = UraniaSunriseSunset.SinAltSun(ddate, hour - 1.0, glong, cglat, sglat) - sinho[j];
			if (ym > 0.0)
			{
				above = 1;
			}
			//
			//  the while loop finds the sin(alt) for sets of three consecutive
			//  hours, and then tests for a single zero crossing in the interval
			//  or for two zero crossings in an interval or for a grazing event
			//  The flags rise and sett are set accordingly
			//
			while (hour < 25 &&(sett == 0 || rise == 0))
			{
				yz = UraniaSunriseSunset.SinAltSun(ddate, hour, glong, cglat, sglat) - sinho[j];
				yp = UraniaSunriseSunset.SinAltSun(ddate, hour + 1.0, glong, cglat, sglat) - sinho[j];
				UraniaSunriseSunset.quad(ym, yz, yp, ref nz, ref z1,ref z2, ref xe, ref ye);
				// case when one event is found in the interval
				if (nz == 1)
				{
					if (ym < 0.0)
					{
						utrise = hour + z1;
						rise = 1;
					}
					else
					{
						utset = hour + z1;
						sett = 1;
					}
				} // end of nz = 1 case
				//
				//   case where two events are found in this interval
				//   (rare but whole reason we are not using simple iteration)
				//
				if (nz == 2)
				{
					if (ye < 0.0)
					{
						utrise = hour + z2;
						utset = hour + z1;
					}
					else
					{
						utrise = hour + z1;
						utset = hour + z2;
					}
					rise = 1;
					sett = 1;
				}
				//
				//   set up the next search interval
				//
				ym = yp;
				hour = hour + 2.0;
				
			} // end of while loop
			//
			// now search has completed, we compile the string to pass back
			// to the user. The string depends on several combinations
			// of the above flag (always above or always below) and the rise
			// and sett flags
			//
			if (rise == 1 || sett == 1)
			{
				if (rise == 1)
				{
					if (EventType > 0)
					{
						OutString = hrsmin(utrise);
					}
				}
				else
				{
					if (EventType > 0)
					{
						OutString = NoEvent;
					}
				}
				if (sett == 1)
				{
					if (EventType < 0)
					{
						OutString = hrsmin(utset);
					}
				}
				else
				{
					if (EventType < 0)
					{
						OutString = NoEvent;
					}
				}
			}
			else
			{
				if (above == 1)
				{
					OutString = AlwaysUp;
				}
				else
				{
					OutString = AlwaysDown;
				}
			}
			return OutString;
		}
		
		public static string sunrise(DateTime ddate, double tz, double glong, double glat)
		{
			//
			//   simple way of calling sunevent() using the Excel date format
			//   returns just the sunrise time or NULL if no event
			//   I used the day(), month() and year() functions in excel to allow
			//   portability to the MAC (different date serial numbers)
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cSunrise);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}
		
		public static string sunset(DateTime ddate, double tz, double glong, double glat)
		{
			//
			//   simple way of calling sunevent() using the Excel date format
			//   returns just the sunset time or ****, ...., ---- as approptiate in a string
			//   I used the day(), month() and year() functions in excel to allow
			//   portability to the MAC (different date serial number base)
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cSunset);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}
		
		public static string CivilTwilightStarts(DateTime ddate, double tz, double glong, double glat)
		{
			//
			//   simple way of calling sunevent() using the Excel date format
			//   returns just the start of civil twilight time or ****, ...., ---- as approptiate
			//   I used the day(), month() and year() functions in excel to allow
			//   portability to the MAC (different date serial numbers)
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cBeginCivilTwilight);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}
		
		public static string CivilTwilightEnds(DateTime ddate, double tz, double glong, double glat)
		{
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cEndCivilTwilight);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}
		
		public static string NauticalTwilightStarts(DateTime ddate, double tz, double glong, double glat)
		{
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cBeginNautTwilight);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}
		
		public static string NauticalTwilightEnds(DateTime ddate, double tz, double glong, double glat)
		{
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cEndNautTwilight);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}
		
		public static string AstroTwilightStarts(DateTime ddate, double tz, double glong, double glat)
		{
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cBeginAstroTwilight);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}
		
		public static string AstroTwilightEnds(DateTime ddate, double tz, double glong, double glat)
		{
			//
			string EventTime;
			string sOut;
			sOut = sunevent(ddate.Year, ddate.Month, ddate.Day, tz, glong, glat, UraniaSunriseSunset.cEndAstroTwilight);
			if (sOut == "....")
			{
				EventTime = "Always Down";
			}
			else if(sOut == "----")
			{
				EventTime = "No event";
			}
			else if(sOut == "****")
			{
				EventTime = "Always Up";
			}
			else
			{
				EventTime = sOut.Substring(0, 2) + "h " + sOut.Substring(2, 2) + "m";
			}
			return EventTime;
		}	
	}
}
